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A Bose-Einstein condensate (BEC) is a quantum fluid that gives rise to interesting shock wave 
nonUnear dynamics. Experiments depict a BEC that exhibits behavior similar to that of a shock 
wave in a compressible gas, eg. traveling fronts with steep gradients. However, the governing 
Gross-Pitaevskii (GP) equation that describes the mean field of a BEC admits no dissipation hence 
classical dissipative shock solutions do not explain the phenomena. Instead, wave dynamics with 
small dispersion is considered and it is shown that this provides a mechanism for the generation of a 
dispersive shock wave (DSW). Computations with the GP equation are compared to experiment with 
excellent agreement. A comparison between a canonical ID dissipative and dispersive shock problem 
shows significant differences in shock structure and shock front speed. Numerical results associated 
with the three dimensional experiment show that three and two dimensional approximations are in 
excellent agreement and one dimensional approximations are in good qualitative agreement. Using 
one dimensional DSW theory it is argued that the experimentally observed blast waves may be 
viewed as dispersive shock waves. 



I. INTRODUCTION 

It is well known that a shock wave in a compressible 
fluid is characterized by a steep jump in gas velocity, 
density, and temperature across which there is a dissipa- 
tion of energy due fundamentally to collisions of particles. 
The aim of this article is to present experimental and nu- 
merical evidence of a different type of shock wave which 
is generated in a quantum fluid that is a Bose-Einstein 
condensate (BEC). In this case, the shock front is dom- 
inated not by dissipation but rather dispersion. Viewed 
locally, these dispersive shock waves (DSWs) with large 
amplitude oscillations and two associated speeds bear 
little resemblance to their classical, dissipative counter- 
parts. However, we demonstrate that a direct compari- 
son is possible when one considers a mean field theory, 
corresponding to the average of a DSW. 

Since extensive theoretical work has been done in the 
field of compressible gas dynamics (cf. Q), it is impor- 
tant to relate this work to the "dispersive gas dynam- 
ics" which BEC embodies. The present work contrasts 
and compares dissipative and dispersive shock waves 
through multidimensional numerical simulation and ana- 
lytical studies. We also provide an explanation of what a 
BEC shock wave is in the context of the well understood 
concept of a classical shock wave in gas dynamics. 

Early experiments studying shock-induced dynamics 
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in BEC were reported in |3| where a slow-light technique 
was used to produce a sharp density depression in a BEC. 
Direct experimental imaging of BEC blast waves has been 
performed in the rotating context and numerical solu- 
tions of the governing Gross-Pitaevskii (GP) equation 
were used to describe the wave dynamics Q ■ Theoretical 
studies of the zero dissipation limit of classicalgas dy- 
namics as applied to BEC was discussed in 0, U where 
it was shown that a shock wave could develop. Subse- 
quently, in the shock wave in the small dispersion 
limit of the one dimensional repulsive GP equation was 
analyzed using the Whitham averaging method and the 
attractive GP equation was analyzed in . 

In the present paper, for the first time a comparison 
between dissipative and dispersive shock waves is carried 
out through a careful investigation of new experiments 
and theory in one, two, and three dimensions. 

The outline of this work is as follows. In section \I~K\ 
we give the relevant dynamical equation, i.e. the GP 
equation, and we put it in non-dimensional form and 
give the associated conservation laws. In section ^ we 
present new experimental results depicting " blast" waves 
in a non-rotating BEC. In section Ulll we show that di- 
rect three dimensional numerical simulations, with radial 
symmetry using the GP equation give excellent agree- 
ment with these experiments. In section Hvl an analysis 
of dissipative and dispersive shock waves in two types of 
one-dimensional systems, the inviscid Burgers' equation 
and the Euler equations, is provided. Two types of lim- 
iting behavior for conservation laws, the dissipative regu- 
larization (small dissipation limit) and dispersive regular- 
ization (small dispersion limit) are considered. We then 
present numerical evidence showing that the three dimen- 
sional and two dimensional calculations agree extremely 
well (less than one percent relative difference). It is also 
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found that the one dimensional approximation of the 3D 
blast wave experiments is in good qualitative agreement. 
Using one dimensional theory, we explain why the exper- 
imentally observed blast waves may be viewed as DSWs. 

A. Gross-Pitaevskii and the Navier-Stokes 
Equations 

An analogy between the classical equations of fluid 
flow, the Navier-Stokes (NS) equations, and the density, 
phase equations for the wave function (order parameter) 
associated with a BEC is well known |^. The crucial 
difference from NS is a dispersive term that replaces the 
dissipative term in classical fluid dynamics. 

The GP equation models the mean-field dynamics of 
the BEC wave function and has been shown to be an 
effective approximation in many situations. Experiments 
in rapidly rotating BECs have provided evidence of dy- 
namics similar to what is often considered to be "blast 
waves" . Moreover, simulations with the Gross-Pitaevskii 
equation were compared with experiment giving support 
to the validity of the GP equation in such extreme cir- 
cumstances Q. 

The dimensional GP equation is 

2 m 

with conservation of particle number 

/ |*pd3a; = l. (1.2) 

The coefficient of nonlinearity, NUq = NAirh'^as/m, is 
characterized by the inter-particle scattering length as 
(here positive representing repulsive particles) and the 
number of condensed atoms N; the other parameters 
are the atomic mass of the species considered (m) and 
Planck's constant divided by 2tt (h). The standard con- 
fining harmonic potential (trap) is given by 

Vo{x, y, ^) = Y {^lix^ + 2/') + ^zZ^) , 

where and uj^ are the radial and axial trap frequen- 
cies respectively. A convenient normalization for our pur- 
poses is to take 

I \ m-^ujl^ ) 

After dropping primes, equation (|l.ip becomes 

= -yV2* + K,^' + l^p^- , (1.3) 

where 



with the normalization of the wavefunction (|L2|) pre- 
served and the trap potential 

Vb(r, z) = i(r^ + a^z^), r^ = a;^ + y^. 

The coefficient = [l^zI^l)^ represents the asymmetry 
in the harmonic trap. In the experiments considered, the 
parameters are: A'^ = 3.5 • 10^ particles, = 5.5 nm 
and TO = 1.45 • 10~^^ kg for the species ^^Rb, [uj^^lOz) = 
27r(8.3, 5.3) Hz, and az ~ 2.45. This normalization shows 
that the dispersion is extremely small, s = 0.012. 

Conservation of "mass" and "momentum" for the GP 
equation l|I.3|) are, 

/ ll'Pdf^^O, 4-/" (**V* - W^-ndf = 0. 

(1-4) 

Since it will be useful for later discussions, we giv e the 
local conservation laws in the ID case (see e.g. [l^) 

Pt + {pu)x = 

£2 (1-5) 

{pu)t + {pv^ + \p^)^ = —{p{\ogp)xx)x - pT4, 

where subscripts denote differentiation. The condensate 
"density" p and "velocity" u are defined by 

* = vpe^^/^ « = 

Equations ljl.5p give an alternative formulation of the 
GP equation in terms of "fluid-like" variables. Since the 

term is obtained from the linear dispersive term in eq. 
(II. 3p . we call this the dispersive term. 

The Navier-Stokes equations for a ID compressible gas 
with density p and velocity u can be written pTl | 

Pt + {pu)x = 

(1.6) 

{pu)t + {pu^ + P)x = ^^^^^ + 

where P is the pressure and F is an external force per unit 
mass. The positive coefficient represents dissipative 
effects due to viscous shear and heat transfer. If the 
pressure law 

P=\P' 

is assumed (for example, a perfect, isentropic gas with 
adiabatic constant 7 = 2 or, equivalently, the shallow 
water equations for height p and velocity u), then the 
NS equations l|I.6|l correspond to the GP conservation 
equations (|I.5|1 when = and F = —Vx- Equations 
(|I.6|1 for the case = and F = are called the Euler 
equations. To compare different types of shock waves, 
we are interested in the dispersive regularization (e^ ~> 
in (|I.5p ') as compared to the dissipative regularization 
(e^ in |jL6l)) of the Euler equations (I4 = F = 0). 
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A remark regarding the form of the above equations. 
Many authors use the velocity form 

Pt + {pu)x = 

Ut + {^U^ + p)^ = —{p{\0gp)a:x)x " 

of equation Ijl.Sp by impUcitly assuming differentiabihty 
and that p ^ 0. As we wiU be interested in weak, hence 
not everywhere smooth, solutions to the dissipative regu- 
larization of the Euler equations, it is necessary to main- 
tain the form of the conservation laws as derived directly 
from the original integral formulation l|I.4(l . It is well 
known that weak solutions to different forms of the same 
conservation law can be quite different. One must also 
be careful when dealing with the vacuum state p = as 
in classical gas dynamics. The momentum equation H1.5() 
takes both of these issues into account. 



II. EXPERIMENT 

In order to investigate the fundamental nature of shock 
waves in a quantum fluid, we have performed new exper- 
iments that involve blast pulses in BECs. In contrast to 
the experiments described in |^, the experiments ana- 
lyzed in this work are all done with non-rotating conden- 
sates. We have succeeded in directly imaging dispersive 
shock waves in these systems, and the particular geom- 
etry of these experiments makes them amenable to the 
theoretical analysis presented in this paper. 

Condensates consisting of approximately 3.5 million 
Rb atoms were prepared in an axisymmetric trap with 
trapping frequencies of {u!±,ujz) = 27r(8.3, 5.3) Hz; uj± is 
the radial frequency and lu^ is the axial frequency. After 
the condensate was formed, a short, tightly focused laser 
beam was pulsed along the z-axis through the center of 
the BEC. The wavelength of the laser was 660 nm, which 
is far red-detuned from the Rb transitions. The pulse 
rapidly pushes atoms from the center of the BEC radi- 
ally outward, leading to the formation of a density ring. 
Before imaging, an anti-trapping technique was used to 
enlarge the features of the blast wave. In brief, a rapid 
expansion of the BEC is created by changing the inter- 
nal state of the atoms such that they are radially expelled 
by the strong magnetic fields forming the trap. Details 
about the anti-trapped expansion are described in [T^ . 
While the anti-trapped expansion changes the scale of 
the features involved, it does not alter the qualitative 
appearance of the shock phenomena, as is confirmed by 
our numerical simulations. 

A sequence of five images taken at the end of experi- 
mental runs with different laser pulse intensities is shown 
in Fig. ^a-e). For this sequence, a 5 ms long pulse was 
sent through the BEC center directly before the start of 
a 50 ms long anti-trapped expansion. The laser waist 
was 13.5 microns. For comparison, the diameter of the 
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FIG. 1: Absorption images of blast pulse experiments with a 
BEC. a-e) Pulse applied before expansion, f) Pulse applied 
during expansion. 



BEC in the radial direction was approximately 65 mi- 
crons. The laser power is given in the images. All im- 
ages were taken at the end of the anti-trapped expansion 
along the z-axis, which is also the direction of the blast 
pulse. For weak blast pulse intensities (Fig. njL,b), essen- 
tially one broad ring of high density is seen, which is due 
to the fact that the laser pulse has pushed atoms radially 
outwards. When the blast pulse intensity is increased, a 
system of many concentric rings appears l^-e) . 

The outcome of a second type of experiment is shown 
in Fig. nf. By pulsing the blast laser during (instead of 
before) the anti-trapped expansion, we can image a situ- 
ation where the compressional ring has not run through 
the condensate yet. For this image, a 5 ms long pulse 
with a power of 1.9 mW and a beam waist of 20 mi- 
crons was used, starting 9.2 ms after the beginning of 
a 55 ms long anti-trapped expansion. In this case, an 
oscillatory wave structure is seen on the outside of the 
compressional ring. The analytical discussion together 
with the numerical studies presented in this paper reveal 
that for both experiments the oscillatory wave structure 
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varying potentials respectively 



FIG. 2: Two examples of experiments with blast pulses in 
slowly rotating BECs. 



is a direct consequence of dispersive shock waves which 
are fundamentally different from classical shock waves. 

Finally, we note that the peculiar wedge shaped ap- 
pearance of the central BEC region in Fig. ^c-e) is due 
to a slight, unavoidable deviation of the laser beam shape 
from cylindrical symmetry. This asymmetry also leads 
to the slightly elliptical appearance of the whole BEC in 
these images. 

By using rotating instead of static condensates, we can 
also observe an intriguing alteration of the blast wave 
pattern. Blast wave images in slowly rotating BECs are 
shown in Fig. [3 Upon slow rotation, dark radially di- 
rected spokes appear in the condensate, cutting through 
the ring shaped pattern familiar from the non-rotating 
case. The number of these spokes increases with increas- 
ing rotation rate, so it is suggestive to attribute them to 
the presence of vortices in the rotating BEC. We specu- 
late that the spokes come about when vortices are present 
in the compressional ring formed by the blast. As this 
ring expands and forms the concentric ring system, the 
density depressions of the vortices are not filled in due 
to the predominantly radial expansion of the compres- 
sional ring. In rapidly rotating BECs, the presence of 
strong Coriolis forces lead to a rather different appear- 
ance of blast waves. The rapidly rotating situation was 
discussed in Q. 



III. SIMULATIONS 



We have performed direct numerical simulations of the 
CP equation (|I.3|I . modeling the two types of experiments 
without rotation explained in the previous section. 

The two experiments depend on when the laser is 
pulsed. We refer to these cases as either in trap (it) 
or out of trap (ot) and model them by the following time 



Vit (r, z,t) 



i(r2 -t- a,z^) + l^e-'-'/''" < t < 6t 
^{-arr^ + a^z^) St < t 



(III.l) 



Vot{r,z,t) = < 



\{r^+a,z^) 



t < 
0<t<U 



^(-arr^ -I- a^z^) + St < t. 



St 



(IIL2) 



The in trap potential (jllLip models a steady state solu- 
tion held in trap while a Gaussian laser pulse is applied 
for the time St. After the laser pulse, a radial anti-trap 
potential is applied [T^ . This models the experiments 
shown in Fig. QJa-e). For the out of trap potential (|III.2|) 
(modeling Fig. QJ), a steady state solution is expanded 
radially. At the time , a Gaussian laser pulse is applied 
with duration St followed by continued radial expansion. 
The Gaussian laser has width d and intensity propor- 
tional to P. The out of trap potential has the effect of 
generating an outward, radial velocity in the BEC before 
the laser pulse is applied whereas the in trap potential 
does not. 

Modeling the non-rotating experiments presented in 
the previous section gives the parameter values: = 
0.48 (9.2 ms), St =^ 0.26 (5 ms), = 0.71 (frequency 
of inverted harmonic potential used for expansion 2tt ■ 7 
Hz), and = 0.57. 

The steady state solution is approximated well by the 
Thomas- Fermi wavefunction 0|, a balance between the 
harmonic trapping potential and nonlinearity. However, 
its use numerically gives rise to unphysical oscillations. 
Therefore, we used an iterative technique similar to the 
technique discussed in to find the precise, smooth 3D 
solution of the CP equation. We provide a brief outline 
of the method here. 

Assume a stationary solution of the form 



*(f,t) 



(in.3) 



where ^ is the normalized condensate chemical potential. 
Inserting the ansatz (|III.3|) into equation (|I.3|I and taking 
its Fourier transform (denoted by ) gives 



(m 



yo0+|0P0=i^(0), 



(ni.4) 



where k is the Fourier wave vector. Equation (|III.4p sug- 
gests the iteration 

7 ^(<^n) 



i£2|fc|2 
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For > 0, the denominator is singular when = 
l^ije^. Also, the normalization is not conserved 

by this scheme. Therefore, wc introduce the modified 
iteration scheme 




(III.5) 

The positive constant c is introduced by adding to 
both sides of equation (|III.4|I . The chemical potential, 
/i, is updated along with the mode (/) by integrating 
equation l|III.4|) . with = (/in+i, multiplied by (0„+i)*. 
This scheme preserves the normalization condition (|I.2|I . 
Given an initial guess 0o ^^nd a large enough value for c 
(we took c = 1.75), we find that the scheme ljIII.5(l con- 
verges to a 3D ground state wavefunction <i>(x) for the 
GP equation. This is a general, fast method for finding 
the 3D BEG ground state for arbitrary potentials Vq. 

The condensate is evolved according to a pseudo- 
spectral Fourier code with a standard 4'^ order Runge- 
Kutta time stepper. Radial spatial derivatives are ap- 
proximated by taking the fast Fourier transform of ^ 
in the radial direction r evenly extended (^'(— r, z,t) = 
5'(r, z,i)) and then multiplying by ik (fc is the Fourier 
wavenumber). The result is then inverted using the 
inverse fast Fourier transform giving a fast, spectrally 
accurate approximation to ^IJj. The zz term is 
approximated similarly. The specific, model equation 
(|I.3|I assumes cylindrical coordinates with axial symme- 

try: V = + r IF ^ 'Ez?- '^^'^ simulations, wc used 
the laser width as an effective fitting parameter. We 
found that by taking the width, d, to be 1.5 times larger 
than its experimental value in the in trap case, excellent 
results were obtained. For the out of trap case, the laser 
width was taken to be half its experimental value, also 
giving excellent results. 

First we consider the in trap case corresponding to the 
potential Ijlll.ip . A plot of the evolution of the conden- 
sate as a function of radial distance is shown in Fig. O 
On the left is the normalized density as given by the 
square modulus of the wavefunction plotted in the 2; = 
plane |5'(r, z = 0,t)p. On the right is the phase gradi- 
ent ^arg^'(r, z = 0,t) or radial velocity in the z = 
plane. The normalized laser width and power used were 
dit = 0.41 and Pn = 0.70 (corresponding to a laser waist 
of 20 ^m and a power of 0.46 mW as in Fig. QJl). 

Fig. El describes the following evolution. The con- 
densate forms a high density ring (i = to 5 ms = 5t 
in non-dimensional form) due to the applied laser pulse. 
When the ring is steep enough (< = 7.5 ms), oscillations 
develop on the inner side of the ring (t = 10 ms). This 
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FIG. 3: BEG density evolved as in the in trap experiment. On 
the left is the density |^'(r, 0, and on the right is the radial 
velocity ^ arg '^{r, 0, t); both are imaged in the z = Q plane. 
The density scale is not constant throughout the sequence; 
in each frame, the density is scaled to its largest value so 
that the dispersing wave is visible. A high density ring forms 
accompanied by oscillations on its inner side, a DSW. This 
DSW expands and propagates outward. Note the time is 
in the center. All figures depicting a time evolution of the 
condensate in this work follow the same form as this figure. 



oscillatory region expands radially due to the inversion 
of the trap potential to an anti-trap potential. The ex- 
pansion continues until time t ~ bb ms when an image is 
taken to compare with experiment (Fig. 

In Fig. 01 we show that the numerical simulation and 
the in trap experiment presented in the previous section 
(Fig. ^) are in good qualitative agreement. The conden- 
sate features at t = 20 ms in Fig. ^expand due to atomic 
repulsion and the anti-trap potential giving the contour 
plot of the density at t = 55 ms (corresponding to the 
end of the experiment) in Fig. 01 The experimental pic- 
ture is reproduced in Fig. 01 left for convenience. The 
simulation used all nominal values for parameters except 
the laser waist {d in non-dimensional units) which was 
taken to be 20 /^m rather than the experimental value of 
13.5 /xm. We speculate that the difference might be due 
to the slight deviation of the experimental beam profile 
from a perfect Gaussian, as indicated by the asymmetry 



FIG. 4: Comparison of the condensate density from the in 
trap experiment (left) and numerical simulation (right) using 
the potential Vu, equation iW.7i . The image from simulation 
is a contour plot of the function J \'i>{r, z,t)\'^ dz, modeling 
the experimental imaging process where the photo was taken 
along the z axis. Both the simulation and experimental pic- 
tures were taken at t = 55 ms. The approximate diameter 
of the condensate cloud from simulation is 850 ^m and from 
experiment is 775 nm. 



of the central regions in Fig. ^ 

Next we consider the out of trap case corresponding 
to the potential Vot pn.2p . A plot of the evolution of 
the condensate in the z = plane as a function of radial 
distance r is shown in Fig. |31 Initially, the condensate 
expands in the radial direction due to the anti-trap poten- 
tial (t = to 9.2 ms or t = in normalized units). A ring 
of high density forms while the laser is on (t = 10 ms), 
similar to the previous in trap case. When the steepness 
of this ring is large enough, oscillations start to develop 
on the inner and outer sides of the ring {t — 10.2 ms). 
These two oscillatory regions expand, quickly overlap- 
ping one another, giving rise to more complicated multi- 
phase type behavior with a propagating wave front that 
continues out radially {t > 10.4 ms). This behavior is 
due to the initial velocity imparted to the condensate by 
the anti-trap potential. The normalized laser width and 
power were dot = 0.21 and Pot = 2.88 corresponding to 
a laser waist of 10 ^m (half the experimental value of 20 
fim) and a power of 1.9 mW. 

In Fig. El wc show that the numerical simulation of the 
out of trap experiment presented in the previous section 
(Fig. nf) show good qualitative agreement. For this com- 
parison, the numerically determined condensate density 
is shown at the time t = 18.0 ms, the result of continued 
expansion from the state at the bottom of Fig. [S] The 
experimental picture was taken at i = 55 ms. 

The rest of this paper is concerned with understanding 
the oscillatory structures that developed in Figures|21and 
|S1 as we will argue that they are dispersive shock waves 
with the oscillations in the latter Figure caused by the 
interaction of two DSWs. 
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FIG. 5: BEG density evolved as in the out of trap experiment. 
On the left is the density and on the right is the radial velocity, 
both imaged in the z = plane. A high density ring forms 
accompanied by oscillations on its inner and outer sides, two 
DSWs. These DSWs quickly interact and propagate outward. 




FIG. 6: Comparison of the condensate density from the out 
of trap experiment (left) and numerical simulation using the 
potential Vot 1III.2I I (right). The image from simulation is 
a contour plot of the function J \'i'{r,z,t)\'^ dz, modeling the 
experimental imaging process. The approximate diameter of 
the condensate cloud from simulation is 116 /im and from 
experiment is 363 pm. 
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IV. CLASSICAL AND DISPERSIVE SHOCK 
WAVES 



In this section, the notions of classical shock waves and 
dispersive shock waves are discussed. We provide a the- 
oretical basis for the experimental and numerical results 
presented in the previous sections. Using the classical gas 
dynamics analogy, it is shown that shocks in a BEC are 
fundamentally different from those in the classical case. 
The analytical methods to understand classical shocks 
and DSWs in one dimension are presented along with 
explicit results for DSW speeds and shock structure. Fi- 
nally, detailed numerical investigations in ID, 2D, and 
3D are presented to show that the qualitative behavior 
of a DSW in 3D is captured by the ID case and that the 
experiments presented earlier do give rise to DSWs. 

We begin by studying the simplest nonlinear dissipa- 
tive and dispersive equations that, under a suitable limit, 
give rise to shock waves, the Burgers and Kortcweg- 
deVries equations. This provides a foundation to under- 
stand the differences between dissipative and dispersive 
shocks. The theory of classical shock waves is well devel- 
oped (cf. P, 0, 01 ) so we will give only a brief synopsis 
suitable for comparison with the much less developed dis- 
persive shock wave case. 



A. Classical Shock Waves, Burgers' Equation 

The classical dissipative shock in one dimension is 
modeled by Burgers' equation 



ut + ilu^)x = e'u^^, (IV.l) 

where is a measure of the dissipation and u represents, 
for example, a density. Equation (|IV.1|I admits traveling 
wave solutions with a hyperbolic tangent profile (see Fig. 
0) 



u{x, t; e) = ^ + ^ tanh { — ^(x — • 



(IV.2) 



The speed of this wave is Vc = 1/2, independent of e. 
In the limit e ^ 0, the (smooth) tanh profile converges 
pointwise to the discontinuous function 



lim u{x, t; e) — u{x, t) — u{x/t) 



The above formula is a mathematical description of a 
classical shock wave. The limiting process e — > 0, e 7^ 0, 
in equation ljIV.l|l is a dissipative regularization of the 
conservation law 




(IV.4) 



Ut + i^u^)x = 0. 



The initial value problem (IVP) for equation pV.4|l 
with the initial data 



u{x, 0) = 




FIG. 7: Dissipative Burgers' equation shock solution <iV.2ll 
with 0, converging to a traveling discontinuity or clas- 

sical shock wave. 



is not well posed because the spatial derivative is un- 
defined at the origin. To see this, note that equation 
(IIV.4P shows that a wave u{x, t) propagates with a speed 
equal to the value of u at that point. Initially, for x < 0, 
the speed is 1 whereas for x > 0, the speed is so w will 
break or become multi- valued at the origin for any t > 0. 
The classical theory of shock waves remedies this prob- 
lem by considering weak solutions and invoking a jump 
or Rankine-Hugoniot condition at a discontinuity which 
relates the value of u on the left (u;) and the right {ur) 
to the speed v of the shock wave. A weak solution u{x, t) 
for the conservation law ljIV.4p satisfies the integral for- 
mulation 

d f'' 

— u(x,t)dx+^(u(b,tf ~u(a,tf) = 0, (IV.6) 

dt J a 

for any a, b such that —00 < a < 6 < 00, thus allowing 
discontinuities in u(x,t) as a function of x. The jump 
condition for Burgers' equation, derived from 



JTVHt as- 
suming a uniformly traveling discontinuity with values 
ui and Ur to the left and right of the discontinuity re- 
spectively, is w = {ui + Ur)/2 or V = ^ for the initial 
data (|IV.5|I . which is exactly the speed of the Burgers 
shock (jIV.2|) . This simple example suggests that finding 
the dissipative regularization of ljIV.4p (equation 
with £ 



HIV.4P (equation IjlV.ip 
0) is equivalent to solving the conservation law 
(|rO|l with the jump condition v = (ui + Ur)/2 at each 
discontinuity. Indeed, this is generally true, assuming the 
entropy condition ui > u,. is satisfied 0| . 

When the entropy condition at a discontinuity is not 
satisfied, u/ < u^, a shock wave solution is not appro- 
priate. The correct choice is a rarefaction wave which is 
continuous for t > 0. Assuming that the solution depends 
on the self-similar variable ^ = x/f, equation 
comes 



'(^-0=0. 



(IV. 7) 



1 X < 
X > 



(IV.5) 



Solutions to ljIV.7p are constants or u{x,t) = x/<, the 
latter corresponding to a rarefaction wave. Then the 
weak solution for initial data u{x, 0) = 0, x < 0, 
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u{x, 0) = 1, a; > is 



u{x, t) — uix/t) 



x/t < V 



J v~ < x/t < v'^ 



(IV.8) 



< x/t 



This rarefaction wave has two associated speeds v~ and 
ti+i = at the interface between the constant left 
state u ~ and the self-similar solution u = x/t and 
= 1 at the interface between the constant right state 
u = 1 and the self-similar solution u = x/t. 

The general case of a system of conservation laws is 
written 

ut + {F{u))^ = 0, 

where u is a vector "density" and F is the vector "flux" . 
The IVP for this n-dimensional system with constant 
step initial data u(x,0) = m;, x < and u{x,Q) = u^, 
X > 0, assuming a dissipative regularization, is called 
a Riemann problem (note that the jump is specified at 
X = 0). The jump condition at a shock with speed v is 



v{Ur - Ui) = F{Ur) - F{ui). 



(IV.9) 



It is well known (see e.g. [l^) that the solution to the 
Riemann problem, assuming certain properties of the flux 
F , is self-similar consisting of n -I- 1 constant states "con- 
nected" by shock waves or rarefaction waves. That is 



u{x/t) = < 



Wo x/t < Vi 

wi{x/t) < x/t < 

Ul < x/t < V2 

Wn{x/t) < x/t < 

U„ V+ < x/t 



(IV.IO) 



where each Uj is constant and Wi{x/t) represents a shock 
or rarefaction wave solution. The Lax entropy condition 
necessary for the existence of the shock wave Wi is 



\i{u. 
Afc(u 
Afc(ui-i) > v., 



1) > 

1) < Vi 



= vt = ^'l > 

and \k{ui) <v, k < i, (IV.ll) 
and Afe(wi) > Vi k > i. 



A shock wave Wi has one speed of propagation so the 2i*'' 
inequality in (|1V.10() is replaced by x = Vit. The Ai in 
(|IV.11|I are the eigenvalues of the matrix with entries at 



fdF,{i 



numbered so that Ai < A2 < ••• < A„. In addition 
to the entropy condition for the i"^ wave Wi to shock, 
the jump condition (|IV.9|) is satisfied for v = Vi 



Ul 



u 1 




FIG. 8: Numerical solutions of equation lIV.12t for initial 
data 1IV.13II and = 0.001 (dashed), = 0.0001 (solid). 
As e decreases, the wavelength of the oscillations decreases. 



Ui-i, and Ur = Ui- Whereas a shock has just one speed, 
associated with every rarefaction wave solution Wi{x/t) 
are two speeds v~ and w,^. 

The established theory of classical shock waves involves 
dissipative regularizations of conservation laws. For ini- 
tial step data, the Riemann problem, there are two types 
of self-similar solutions of interest, a shock wave and a 
rarefaction wave. In the next section, we study what 
happens in the region nearby breaking when a disper- 
sive, rather than dissipative, term is used to regularize 
the conservation law. It will be shown that self-similarity 
plays a crucial role and that a dispersive shock wave cor- 
responds, in some sense, to a simple rarefaction wave 
solution of a system of conservation laws. 



B. Dispersive Shock Waves, Korteweg-deVries 
Equation 

As a simple model of dispersive shock waves (DSWs) 
(e.g. in plasmas we consider the Korteweg-deVries 

(KdV) equation 



(IV.12) 



for £^ ^ 1. We investigate the behavior of the solution 
to the initial value problem (IVP) 



u(a;,0;e) 



1 a; < 
a; > 



(IV.13) 



as e — > 0. This is a Riemann problem in the context of 
a dispersive regularization of the conservation law ljIV.4|) 
with no inherent dissipation. 

Figure [S] depicts two numerical solutions to the IVP 
(IIV.12|) and (|IV.13p for smaU e^. OsciUations develop 
about the initial discontinuity with wavelength propor- 



tional to e. Then, as e 



0, an infinite number of 
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oscillations develop. To understand this behavior, one 
can employ Whitham's method ^3 which is an exten- 
sion of the Krylov/Bogoliubov method of averaging for 
ODEs to PDEs. The essence of the technique is to as- 
sume that the KdV equation l|IV.12|) has a uniformly 
traveling wave solution and average the PDE's conserva- 
tion laws over fast oscillations allowing certain parame- 
ters (such as amplitude, wavelength, and speed) to vary 
slowly in time and space. Gurevich and Pitaevskii |23| 
applied Whitham's method to the IVP considered here 
with boundary matching where one derives boundary 
conditions for the oscillatory region based on continu- 
ity of the averaged flow. We will follow Gurevich and 
Pitaevskii's work, applying the method of initial data 
regularization presented in |2ll | rather than boundary 
matching to asymptotically solve the initial value prob- 
lem (|IV.12|) and (|IV.13|I in the hmit ^ 0. This hmit, 
denoted u, is a weak limit where one averages over the os- 
cillations and is different from the Burgers' shock strong 
limit depicted in Fig. [7| The Whitham method describes 
the asymptotic (large t) behavior of the slowly modulated 
oscillatory region seen in Fig. |S1 which is a dispersive 
shock wave, by enabling an explicit calculation of the 
weak limit u. 

Just as a discontinuity represents an idealized dissipa- 
tive shock wave in a compressible fluid, the weak limit u 
represents an idealized dispersive shock wave. Any com- 
pressible gas will have a small but non-zero amount of 
dissipation. Since a strong limit exists, the transition 
from small to zero dissipation is smooth. In the case 
of a DSW, where a weak limit prevails, the transition 
from small to zero dispersion is accompanied by an in- 
finite number of oscillations. Thus, any physical DSW 
with small but non-zero dispersion will consist of a finite 
number of oscillations. However, the weak limit u pro- 
vides an understanding of the physical DSW structure 
and its associated speeds. As we will show, it also en- 
ables clear comparisons between classical and dispersive 
shock waves. 

The first step in the Whitham averaging method is to 
obtain a quasi-stationary periodic solution. Assuming 
the traveling wave ansatz, u{x,t;e) = <t'{())^ 9 ~ (x — 
Vt)/e, equation HIV.12|) reduces to 

-V(j3' + + (j)'" = 0. 

Integrating this ODE twice we obtain 

1 , , 



A'^2 



+ A(t) + B) = \P{(t)) , 



with A and B arbitrary integration constants. Solutions 
to equations of this form, when P{(t>) is a cubic or quartic 
polynomial, arc elliptic functions. We write the polyno- 
mial P in terms of its roots 

P(0) = (Ai - (/))(A2 - 0)(A3 - 0), Ai < A2 < A3. 

For convenience, we make the following linear transfor- 




FIG. 9: Elliptic function solutions to the KdV equation for 
several choices of the parameters {vi}. The solution converges 
to a constant as m ^ and to the soHton sech profile as 
m ^ 1. 



mat ion 

ri + A2), r2 = i(Ai + A3), rs = ^(Az + A3), 

''1 <r2< r-i. 

Then the elliptic function solution can be written as (see 



(t){e) = ri + ra - r3 + 2(r3 - ri) dn 

f2~Tl „ X ^ Vt 



/ r3 - ri 



V 



rs - ri 



6 ' y 

^{n + r2 + 7-3). 

(IV. 14) 



This is an exact solution to equation (|IV.12p with three 
free parameters {ri} related to the amplitude: max(0) — 
min((?!)) = 2(r2 — ri), speed V, and wavelength 



L = 2K(m) 



where K(m) is the complete elliptic integral of the first 
kind. Note that L is obtained from the periodicity of the 

dn fimction \\Y AA} i.e. [^J~^^^0] = 2K{m) where [•] 
is the period of the argument. The parameter m is the 
modulus of the elliptic function. See Fig. I^for a plot of 4> 
for various values of m. There are two limiting behaviors 
dn(?/;0) = 1 and dn(?/; 1) = sech(y), the solitary wave 
solution. 

The basic idea behind Whitham theory is in the pro- 
cess of averaging over "fast" oscillations. This yields 
the behavior of the weak limit. It, of equation ljIV.12p 
as e — !■ 0. Since e is assumed to be much smaller than 
1 in cq. ljIV.12p . the phase 9 = {x — Vt)/£ is a fast 
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variable. We assume that modulations of this periodic 
solution take place on the scale of the "slow" variables x 
and t. Then the average of is 



9,x,t) de 



ri{x,t)+r2{x,t)-r3{x,t)+ (IV.15) 
E[m{x,t)] 



2[r3ix,t)~ri{x,t)] 



K[m{x,t)] 



where E{m) is the complete elliptic integral of the second 
kind. 

The next step is to write down the first three conser- 
vation equations for the KdV equation |23| 







(l"^) + (i"^ + ^"^u^xx - ^£^ul^ = (IV.16) 



2e^uul 



We require three equations because there are three pa- 
rameters {ri\ that arc allowed to slowly vary in time and 
space. Now we insert the periodic elliptic function solu- 
tion (j) into equations l)IV.16(l and average the equations 
over the period L to find 

(l^ - ^) ^ + (i^ - 400f + 30|) ^ = 0. 

Note that (px = (t^e/e- 

Assuming that the parameters {n} depend on the slow 
variables x and the above equations can be transformed 
to Riemann invariant form |l9ll20| 

UT ' UT 

^+f.(r-i,r2,r3)^ =0, i = 1,2,3. (IV.17a) 

The variables {r.i} are the Riemann invariants for the 
hyperbolic system ljIV.17ap with the velocities 

1, .2, . K{m) 

6 6 K[m) — h[m) 

1, ^2, , (l-m)A'(m) 

3 3 L[m) — [1 — ■m)K(m) 

1, .2, Jl-m)K{m) 

^'3 = +r2 +r3) - -(rs -ri)> 



E{m) 



(IV. 17b) 



We wish to solve equations ljIV.17(l subject to the ini- 
tial data l|IV.13|l . This is accomplished by the method 
of initial data rcgularization, presented in [2l| . Although 



1 = 









■■3 












u 











FIG. 10; Initial data regularization for the KdV dispersive 
Riemann problem. The dashed line represents the initial data 
l|IV.13|l for u. The solid lines represent the initial data for the 
Riemann invariants ri, r2, and ra that regularize the initial 
data for u. This initial data for the Riemann invariants sat- 
isfies the three properties of characterization, non-decreasing, 
and ordering IIIV.lQI l so a rarefaction wave solution exists for 
all time (see Fig. Illjl . 



the background to this method involves a detailed under- 
standing of inverse spectral theory and Riemann surface 
theory, the method itself is straightforward. Any solution 



to equation ljIV.4() with decreasing initial data will even- 
tually break. The dissipativc rcgularization handles this 
by introducing jump conditions across the shock relating 
its speed to its values before and after the discontinuity. 
The dispersive regularization employs the higher order 
hyperbolic system ljIV.17p with initial data that char- 
acterizes the initial data for u, is non-decreasing, and 
satisfies a separability condition. The initial data 



ri(x,0) = 0, r2(a;,0) 



X < 



1 x>0 ' '^3(^,0) = 1, 

(IV.18) 

shown in Fig. UHl regularizes the IVP (|IV.12p and HIV.13|) 

because of the following properties 



0(x, 0) = u{x, 0; e) 
fe(^-0)>0 



(characterization) , 
(non-decreasing), 



maxri(a;,0) < minri_|-i(a;, 0) (separability). 

(IV. 19) 

Characterization amounts to verifying that the initial 
data for the full problem ljIV.13() is equivalent to the 
initial data for the averaged problem cj); the same as- 
sumption is made in the boundary matching method 
|20l | . The non-decreasing and separability of the r.i ensure 
that a global, continuous (non-breaking) solution to the 
Whitham equations IjlV.lZp exists for all time [l5ll2^. 
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FIG. 11: Dissipative (dashed) and dispersive (solid) regular- 
izations of the conservation law <1V.4|| . The dissipative case 
corresponds to a traveling discontinuity with speed i satis- 
fying the jump condition lIV.9t . The dispersive case is a 
rarefaction wave solution to the Whitham equations <iV.17ll 
with two associated speeds iij = 2/3 and V2 = — 1. This rar- 
efaction wave modulates the periodic solution (llV.14t giving 
a DSW (see Fig. [Eland eq. HIV.23|l l. 



The system (jlV.lTII with initial data (jlV.lSp has an 
exact rarefaction solution in the form of a self-similar 
simple wave with ri = 0, = 1, and r2{x,t) = r2{£,)i 

= x/t. The remaining nontrivial equation in (jIV.17a(l 
takes the form 

{V2 - = 0, 
which is satisfied when the implicit relation U2 = ^ or 

[l-r2(CM[r2(0] 



i[l + ^2(C)]-^2(C)- 



E[r2m-[^-r2mK[r2m 

(IV.20) 

is satisfied. The above is one equation for one unknown, 
r2(C), which is solved by a standard root finding method 
for each ^ (see Fig. ITT)) . 

The rarefaction wave has two associated speeds and 
V2 which are determined from the Whitham equations 
(I1V.17II . Ahead of the moving fronts, the arc constant. 
Since 



dr2 
~dt 



0, when 



dx 



V2, 



from equations ljIV.17b|l , the speeds are given by the lim- 
its 

v+= lim V2{0,r2,l)^l, (IV.21) 

7-2 ^ 1 ~ O 

v~ = lim W2(0,r2,l) = -1. (IV.22) 

The dispersive Riemann problem, equation ljIV.12|) 
with initial data ljIV.13p or, cquivalcntly, equations 



Classical 




1/2 2/3 



FIG. 12: Comparison of a dispersive shock (eq. I|IV.23^ with 
£^ = 0.0001), its average (f) (eq. 11 V.15L dashed), and a clas- 
sical, dissipative shock (eq. I)IV.3^ . Burgers' solution) plotted 
at time t = 1. The DSW front moves faster than its classical 
counterpart. The average (p looks very similar to the classical 
shock, a steep front connected to a constant in the rear, ex- 
cept that the speed of the front is different and the function 
is continuous. 



(|IV.14|) and HIV.17|) with initial data HIV.18|) . has the 

asymptotic (i » 1 and ^ 1) DSW solution 

u{x,t;e) « r2{x/t) - 1 + 2dn2 ( ^— m = r2{x/t) 



V{x/t) = Ul + r2{x/t)). 



(IV.23) 



The function 7'2 (x/t) is the rarefaction wave solution sat- 
isfying equation (|IV.20|1 . This DSW solution, its average 
(eq. (jlV.lSp ). and the Burgers type classical shock solu- 
tion ljIV.3|) are sliown in Fig. ^1 The DSW averaging 
process produces a shock front that resembles the classi- 
cal shock front but, however, has a different speed and the 
DSW front is continuous. The front speed of the DSW, 
|, is the phase speed of the classical soliton solution to 
KdV which fixes an amplitude of 2 

M(x,t) = 2sech2(-^(x-|i)). 

One can think of a DSW as a slowly modulated train of 
solitons decaying, in a self-similar fasliion, to a constant. 
The DSW is based on the rarefaction solution (|IV.20|) so 
it has two associated speeds, the trailing edge V2 (|IV.22() 
and the leading edge V2 HIV.21|) . Even though the DSW 
is non-zero as a; ^ —00, the oscillations remain in a finite 
region of space describing the expanding behavior of a 
steep gradient in this dispersive system. 

In Fig. 1131 we show the numerical solution to the KdV 
equation with the step initial data (|IV.13p for = 0.001. 
The wavelength of oscillation, leading edge amplitude, 
and speed of the asymptotic solution agree well with the 
numerical result. The position of the leading edges differ 
slightly because the asymptotic solution is valid for i ^ 1 
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KdV Numerics 



1.5 




Asymptotic Solution 




FIG. 13: Numerical solution of the KdV equation with an 
initial step (top) and the asymptotic DSW solution 11 V. 2311 
(bottom) for = 0.001. 



and it takes the numerical solution some time to reach 
this stage. Note that the speed of the leading edge in the 
numerical solution, averaged from t = 5 to i = 7, is 0.660 
which is approximately |, the analytical result (|IV.21|) . 

If the initial data for the KdV equation ljIV.12p is non- 
decreasing then no breaking occurs and a global solution 
exists for the zero dispersion limit (e^ 0). For the case 



0; e) 



a; < 

1 a; > 



(IV.24) 



the solution is a rarefaction wave, a weak solution sat- 
isfying the conservation law 



l|IV.4p . This is the Burgers 
rarefaction wave HIV.8() (see Fig. I14|) . 

Whitham averaging provides an effective way to define 
the DSW shock speeds and derive the asymptotic oscilla- 
tory structure of a DSW along with its leading amplitude. 
We note that in any experiment e will be finite hence os- 
cillations will exist. The averaged solution is useful when 
comparing with gas dynamics since we can evaluate the 
jump in density across a DSW region and determine the 
velocity of a DSW shock front. 

The long time asymptotic behavior of the KdV equa- 
tion has been analyzed in |25| . In general, an arbitrary 





0.05 




0.005 


,.,.,p2_ 






0.5 



-0.5 




-4 



FIG. 14: Numerical solution of the KdV equation HIV. 1211 
with the step initial data l|IV.24^ for different values of e and 
the zero dissipation/dispersion limit e = 0, the rarefaction 
wave l)IV.8|l . The plot corresponds to t — 1. 



initial condition will evolve into a dispersive tail |25l | , a set 
of sohtons [i^ili^l, and a " coUisionless shock" region [25^ . 
In a sense, the asymptotic solution in Fig. ^Icontains all 
of these regions. The very front of the oscillations is the 
coUisionless shock region over which a constant connects 
to a train of sech^ solitons eventually leading to small, 
linear oscillations at the tail. 

Thus we have described how to study a dispersive 
shock wave associated with KdV in the context of 
Whitham theory. A DSW can arise in the dispersive reg- 
ularization of a conservation law just as a classical shock 
can arise in the dissipative rcgularization of a conserva- 
tion law. The key difference is that a weak limit where 
one averages over the oscillations is required in the dis- 
persive case. This method gives useful results such as the 
asymptotic modulated oscillatory profile, the wavelength 
of oscillation, the leading amplitude, and the speeds of a 
dispersive shock. On a large scale, once the limiting pro- 
cess has been accomplished, the DSW and classical shock 
look similar, i.e. constants connected by sharp gradients. 
However, the shock speeds are different. 

On a small mathematical note, Lax and Levermore 
used the inverse scattering transform to take the limit 
e in the KdV equation ljIV.12p for a broad class of 
initial data. They showed that the limit, m, is a weak 
limit in the sense that 



lim 



u{x, t; e)f(x) dx 



u{x, t)f{x) dx 



for all smooth, compactly defined functions f{x). This 
type of limiting procedure is required because the solu- 
tion develops an infinite number of oscillations. Lax and 
Levermore also showed that, in a region of breaking, the 
weak limit u can be calculated explicitly by using the 
Whitham averaging method thus giving the method a 
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strong mathematical footing. 

Throughout the rest of this paper, we will use this dis- 
sipative/dispersive analogy with the Burgers' and KdV 
equations to motivate our discussion of the more com- 
plicated problem involving the dissipative and dispersive 
regularizations arising in the context of BEC and gas dy- 



C. Dissipative Regularization of the Euler 
Equations 

As mentioned in the introduction, the compressible 
equations of gas dynamics without dissipation are the 
same as the local conservation equations for a BEC ljl.5|) 
with e = 0. Let us consider the Ricmann problem for the 
dissipative regularization of the Euler equations in one 
dimension with step initial data 

Pt + {pu)x = 
{pu)t + {pu^ + = 



p(x,0) 



1 .T > 



u(x, 0) 



Uq .t < 
X > 



(IV.25) 



This is a general step initial value problem with two pa- 
rameters po and uq. Note if we make the transformation 



p = prP , u = \fp^U -f Mr , 

t = Prt , X = \f(h-{i — Urt), Pr ^ 0, 



(IV.26) 



then we convert the initial conditions in l|IV.25() to the 
general step initial conditions 



p(i,0) 
u{x, 0) 



Pi = PrPO X <0 
Pr X > 

Ul ~ y/p^UQ + Ur X < 
Ur i > 



The unique weak solution to l|IV.25p consists of, in 
general, three constant states connected to one another 
via two wave solutions: shocks or centered rarefaction 
waves (this is a special case of (jIV.10|) for n ~ 2). We will 
first study the canonical, right-going shock case which 
consists of two constant states connected by a traveling 
discontinuity. 

The jump conditions l|IV.9p for the IVP 1IV.25|I are, 
assuming a traveling wave shock. 



Vipo - 1) = PoUq 
V{pqUq) = poul + \pI 



(IV.27) 



A physically realizable shock solution must satisfy an 
entropy condition, the statement that across any shock 
wave, there must be a corresponding increase in entropy. 
For the Ricmann problem ljIV.25|) , the entropy condition 



for a right going shock is simply po > 1 ^3 • The entropy 
condition and the jump conditions (|1V.27|I determine a 
classical shock wave uniquely. Solving for uq and V in 
(|IV.27|) gives 



"0 = ±(Po - 1) 



(IV.28) 



with the corresponding shock speed 

For an entropy satisfying right-going shock, we take the 
plus sign in ljIV.28|) . With this specific choice for mq, the 
Riemann problem (|IV.25p has the unique weak solution, 
parameterized by po. 



a shock moving with speed 



u{x/t) 



Uq x/t <V 

x/t>V 



V = poi 



I pq - 1/Po 
2(po - 1) ■ 



Note that for a right-going shock to exist, there must be 
a non-zero density p on the right. Otherwise, the solution 
is purely a rarefaction wave |l7l | . 

When Po < 1 ; a pure rarefaction solution exists for the 
specific choice 

Mo 1). 

This choice is a consistency condition for the existence 
of a continuous rarefaction wave connecting the left con- 
stant state (potPqUo) and the right state (1,0) The 
rarefaction solution is given by (also sec Fig. I20() 



pix/t) 



u{x/t) = 



Po x/t<{?,^-2) 
i(2 + x/t)2 (3Vp^-2)<x/t<l 



1 < x/t 



2^ - 2 x/t < (3VaI - 2) 
i(-2 + 2a;/i)2 (3^ - 2) < a;/t < 1 . 

1 < x/t 

(IV.29) 







By manipulation, the Euler equations ljIV.25p can be 
written in the Ricmann invariant form 



dt 



dr- 1 ,„ , dr- 

^ + 1(3-++-) — 



dx 
W-2VP 



(IV.30) 



These equations yield a general solution of the form 
(jlV.lOp . In the next section, we will compare the so- 
lution of equations (|IV.30|) to the zero dispersion limit of 
the GP equation. 
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The dissipative regularization of the Euler equations- 
the jump conditions HlV.27(l "gives a criterion for deter- 
mining the speed of a classical shock wave given the ini- 
tial jump in density. In the next section, we will discuss 
DSWs in BEC and show that rarefaction waves play a 
crucial role in their understanding. 



Dispersive Regularization of the Euler 
Equations, BEC 



D. 



Consider equations ljl.5|) . the local conservation equa- 
tions for a BEC, with 0. This is the zero dispersion 
limit of the CP equation, which wc consider as the dis- 
persive regularization of the Euler equations considered 
in the previous section. Assuming free expansion of the 
condensate or zero potential, = 0, the CP equation is 
equivalent to the 3D, defocusing nonlinear Schrodinger 
equation (NLS). Later we show that the ID NLS equa- 
tion gives good qualitative agreement with the three di- 
mensional problem in the so-called "blast wave" regime. 
In this regard, the Whitham averaging method is a useful 
device to analyze DSWs in a BEC. It has been used to 
analyze the ID defocusing NLS equation in the case of 
an initial jump in density and velocity using boundary 
matching j^, [sO] ■ We apply Whitham averaging to 
the ID NLS equation using initial data regularization as 
presented in |2l|, |31| to derive the shock structure of the 
canonical ID BEC DSW along with its associated speeds. 
This technique is equivalent to the boundary matching 
method for a single dispersive shock wave. It also allows 
for generalizations to more complicated, multi-phase type 
interactions which we will study in the future. 

For small dispersion, oscillations begin to develop in 
breaking regions. As in the KdV case, a strong limit 
does not exist; hence we are lead to consider a weak limit 
where one averages over the oscillations. Consider the 
dispersive Riemann problem for BEC in one dimension 
without an external potential V, which models a freely 
expanding condensate 



Pt + {pu)x = 



{pu)t + {pu^ + -p^)x 



p(a;,0;e) 



4 

(IV.31) 

Recall that p ~ j^Pp represents the condensate density 
and u = £(arg is the condensate flow velocity. These 
initial conditions are general step-like data (see (|IV.26|l '). 

To apply Whitham theory, wc require a periodic solu- 
tion to equation (|IV.3ip . the local conservation equations 
for the ID NLS equation 

1^2, 



I^P^- = 0. (IV.32) 
Assume a traveling wave solution of the form 

= ^(6l)e"^('')/^ e={x-Vt)/e. 



Inserting this ansatz into equation ljIV.32p , and equating 
real and imaginary parts gives 



-VA' + (^'A! \(i>"A = 
V^A + iA" - -A^ = ^. 



(IV.33) 



Integrating the first equation and solving for (/)', wc find 

2ci 



where C\ is a constant of integration. Inserting this result 
into the second equation in ljIV.33p and simplifying gives 

A" -f V'^A - —1 - 2^3 = 0. 
A3 

Integrating this equation gives 

+ 1/2^2 + M _ ^4 ^ C2 = 0, (IV.34) 

where C2 is a second constant of integration. To obtain 
the elliptic function solution, let p ~ A^] then equation 
(IIV.34P becomes 



p'2=4(p3_l/V-2c2p-4c?) 
^ A{p - \i){p - \2){p - 



(IV.35) 



< Ai < A2 < A3. 
The periodic solution to the above equation is 

p{x, i) = A3 - (A3 - Ai)dn2(VA3-Ai ^''^^^^ m) 



u{x,t) = (j)'{9) = V - 



2ci 

p{x,t)' 



V = VAi + A2 + A3 



•^2 - Ai 2 1 \ \ \ 

m^- — , ci = JA1A2A3. 

A3 — Ai 



(IV.36) 



Similar to the KdV equation, this solution has three in- 
dependent constants of integration Ai, i = 1, 2, 3, the 
roots of the cubic polynomial in HIV.35|) . However, by 
the invariance of the NLS equation with respect to the 
"Galilean boost" 

*(a;, t) ^ e-'^("-5^*'*(x - Vt, t), 

the phase speed V is another arbitrary constant. Then 
the periodic solution to equation (|IV.31|) is 

p{x, t- e) = ^{9) = A3 - (A3 - Ai)dn2(VA3 - XiQ- m) 



u{x, t; e) = ^{9) — V — a 



V A1A2A3 

m 



A2 - Ai 

m=- — , cr = ±1, 

A3 — Ai 



x-Vt 



(IV.37) 
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with four arbitrary parameters V and A^, i = 1,2, 3. The 
sign of the constant of integration ci in (|1V.36|I is not 
determined. Either sign a = ±1 gives a vahd periodic 
solution to the NLS equation (|IV.32p . This will be im- 
portant later in our analysis of DSWs with points of zero 
density. 

Using Whitham's averaging method, we will derive 
equations describing slow modulations of the four param- 
eters in the periodic solution l)IV.37p . Anticipating the 
form of the Whitham equations, we will express the four 
aforementioned parameters in terms of four parameters 



V = i(ri +7-2 +r3 +r4), 
A2 = + ^2 - + 

ri < 7'2 < rs < r4. 



(IV.38) 



Similar to the KdV system, the parameters correspond 
to Riemann invariants for the NLS Whitham equations, 
which we derive now. 

The period of the dn function in ljIV.37|) is 
[VA3 — XiO] = 2K{m) where [•] denotes the period of 
the argument. Then the wavelength of oscillation for the 
periodic solution pV.37(l is 



L 



2K{m) 
\/A3 — Ai 



(IV.39) 



The average of the periodic solution ljIV.37p over the 
wavelength L is 



ij{x,t) 



ipi'{x, t) 



v(x, t) 



A (X A^^^™) 



Vip - cr VAiA^, 



L 



iy{9) d9 



A, 



F{x, 1 - m)[E{m)/K{m) - 1]} 



X = sm 



A3 — Ai 



A3 



(IV.40) 



The functions F{x, 1 — m) and E{x, 1 — m) are incomplete 
ellip tic integrals of the first and second kind respectively 

m. 



Four conservation laws for the NLS equation are 0, 



Pt + [pu)x = 



/ 1 

{pu)t + [pu^ + -p^ - —p{\ogp) 







Ap / 1 \ 4p / X 

e'pxxx - Is'P^ + 1^34 - 5eppx - 3epxu' + 

_ 5 ^4 PxxxPx ^ 3 ^2 Pxx , 
xxxx ~r 

p p 



3epuux^^ + (^^e'^P. 



g c ^2 2 yP'-cx 2 rxx gC iJ,j.-r 



-Qe'^uUxPx + 2p-^ + 7p^w^ 



be^puuxx - ^s^pul + 2pu^ + 2p^u^\ = 

/ X 

(IV.41) 



To obtain the Whitham equations for the {ri}, we insert 
the periodic solution ljIV.37p into the conservation laws 
PV.41|1 and average over the fast variable 9 to find 



i3i 

4^- 



..... .. .. .i- 



8 ^2 



(IV.42) 



Assuming that the four parameters vary on the slow 
length and time scales x and t, the Whitham equations 
are obtained Ha, 113 



^+«,(ri,r2,r3,r4)^ =0, z = 1,2,3,4. (IV.43a) 



The Vi are expressions involving complete first and sec- 
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ond elliptic integrals 



vi = V ■ 



V2 = V + 



ri) 



1 



V4 = V + 



^(^•4-^3) 1 

{rj - r3){r2 - ri) 
{ri - r2){rs, - ri)' 



(^4 

v ^ 


— ro)E(m) 


-1 




— r\)K(m) 




V i 




-1 




— r2)K{m)_ 




(r-4 


— r2)E{m) 


-1 


(^3 


— r2)K{m) 






— ri)E{m) 


-1 


(^4 


— ri)K{'m) 





(IV.43b) 



Thus (the complicated looking) equations (|IV.42|I reduce 
to the simple system of first order PDE HIV.43|I . This 
reduction takes advantage of certain integrable properties 
of the NLS equation which we will not discuss here. For 
further details, the reader is referred to |3^ . 

As in the KdV case, it turns out that the right going 
dispersive shock wave is characterized by a self-similar, 
simple rarefaction wave in the associated Whitham equa- 
tions (|IV.43|) . There are two free parameters in the IVP 
(|1V.31|I . the initial velocity uq and the initial density pq. 
As we are seeking a simple wave solution, wc require one 
of the Riemann invariants from the Euler system (jIV.30|l . 
r_ in this case, to be constant initially. For the initial 
data in (|1V.31() , this corresponds to the specific choice for 
the initial velocity uq — 2(y/po — 1) and the initial data 



r+ix,0) 



Po 
2 



.T < 

a; > 



r_(2;,0) = -2. (IV.44) 



Equations pV.30|l accurately describe a regular solu- 
tion, hence a dispersive or dissipative regularization of 
the Euler equations whenever both and r_ are non- 
decreasing |2l|. A weak limit is not required because 
a rarefaction type solution exists for all time. In other 
words, the dissipative and dispersive regularizations are 
equivalent when no shocks develop. 

On the other hand, when the initial data for the Rie- 
mann invariants are decreasing, as in the case of l|IV.44|) 
with Pq > 1, Si shock wave will develop. 

We now use the technique of initial data regularization 
[2l] | where the ri are chosen initially (see Fig. I15|) so that 



'(/'(x, 0) = p(a;, 0; e), iy{x,0) = u{x,0; e) (characterization), 
^(a^jO)>0 (non-decreasing). 



maxri(a;,0) < min 7-^+1 (x, 0) 



(separability). 
(IV.45) 

Recall that the characterization property implies that the 
initial data for the full problem and the averaged problem 



FIG. 15: Regularized initial data for the dispersive Riemann 
problem II1V.3H . The variables r+ and r_ are the Riemann 
invariants for the Euler equations satisfying <IV.30ll . The 
{ri} axe the Riemann invariants for the Whitham equations 
IIIV.43I 1 describing the modulation of a periodic wave (IIV.37I 1. 
The {r,} satisfy the properties of being non-decreasing and 
separable <IV.45ll . so a continuous rarefaction solution exists 
for all time (see Fig. Ilfcil . 



are the same. For this to be so, we take 



ri(x,0) = -2, r2(x,0) = 2, ri{x,0) = ~ 2, 

2 a; < 
4Vpo - 2 a; > 

(IV.46) 



r3(x,0) 



We also require a spatial dependence of the sign a in 
(|IV.40|) when po > i 



(j{x, 0) = sgn(a;) 



-1 a; < 
1 a; > 



(IV.47) 



When 1 < Po < 4, cr = 1. The non-decreasing and 
separability properties guarantee that a continuous solu- 
tion for the hyperbolic system (|IV.43|I exists for all time 
[Him. Thus we have regularized the "shock" initial 
data. 

Using the initial data regularization shown in Fig. 1151 
we solve equations (|IV.43|) for a self-similar (^ = x/t), 
simple rarefaction wave. Assuming ri = —2, r2 = 2, 
r3(a;,t) = T'alC), and r4 — '^^fpa ~ 2, we find that all the 
Whitham equations are satisfied if 



(^^3 - iVz = 0. 



This equation is satisfied for non-constant r3 when vj, 
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FIG. 16: Solution to classical (dashed) and dispersive (solid) 
Riemann problems for the Euler equations with initial data 
shown in Fig. 1151 The classical solution consists of three 
constant states connected by a rarefaction wave and a shock 
wave. The dispersive regularization involves a pure rarefac- 
tion solution of the Whitham equations <lV.43t where only 
ra varies according to <IV.48ll . This solution modulates the 
periodic solution <1V.37II giving a DSW (see Fig. I17II . 



^, or explicitly 




2- 



(r3(0-2)A'(m(0) 



(VA^-l)(r3(^) + 2)- 



(IV.. 



The above is one nonlinear equation for the unknown 
'^3(0- We use a standard root finding method to solve 
this system. The rarefaction solution for is shown in 
Fig. ^] along with the dissipative regularization of the 
Euler equations l|IV.30(l for the same initial data. Re- 
call that the general solution for the dissipative Riemann 
problem has the form (|IV.10|) for n = 2. For the initial 
data in l|IV.44|) . we find that this general solution consists 
of three constant states, the first two connected by a rar- 
efaction wave, the latter two connected by a dissipative 
shock wave (see the dashed curve in Fig. I16|l . 

Inserting the self-similar solution of the Whitham 
equations l|IV.48|) into the original periodic solution 
(|1V.37|I for the case 1 < po < 4 (recall cr = 1), gives 
a slowly modulated periodic wave (see Fig. [T7| . a DSW 
in BEG. 

The speeds of the trailing and leading edges of the 
DSW are found in the same manner as they were in the 



FIG. 17: DSW (GP, osciUatory) and classical shock (Euler, 
discontinuity). The dashed curve is the averaged DSW. The 
upper plot is the density p, the middle plot is the velocity 
u, and the lower plot is the momentum pu for the parame- 
ters po = 2.5, f = 1, and = 0.002. A DSW in BEG is an 
expanding region with the constant trailing and leading edge 



speeds — ^/p^ and 



8po-8^+l 



respectively. When 



■-'3 - 2^-1 

1 < Po < 4, the DSW looks similar to the one pictured here. 
For Po > 4, the situation is different see Fig. 1211 The maxi- 
mums and minimums of the density and velocity arc marked 
for comparison with the NLS gray soliton solution <IV.50t . 



KdV case, using equations (|IV.43bp 

Wg- = lim vz{~2,2,rz,4.^-2) 
v+ = lim W3(-2,2,r3,4V;5^^-2) 
8po-8Vp^+l 



P^^, (IV.49) 



2VP^-1 

Because the trailing edge speed is the speed of sound 
(linear disturbances of the Euler equations (|IV.25(l ) in gas 
dynamics (see e.g. jl7|). it might appear that the trailing 
edge of a DSW is not affected by the nonlinearities in 
the problem. However a closer examination reveals that 
the trailing edge speed is exactly the phase speed of a 
gray soliton with the trailing dip shape given in Fig. 1171 
A general gray soliton solution to the ID NLS equation 
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(II V. 3211 is 

^B(3^ sech^ (61) 



-W 



S2tanh^(6') + ^2 
e = Bf3[x - [Eli - W)t - xo]/e. 

(IV.50) 

where all parameters are real and /i^ + /J^ = 1 (see Fig. 
I18|l . The speed of the gray soliton is 



(IV.51) 



It's minimum density and velocity occur when x = {Bj.i- 
W)t + xq giving the maximum and minimum values 



Pmax = maxl^p = 
Umax = maxe(arg\E')2; = 

Urnin = min£(arg*)2; = B[p 



Pmin = minl^-l^ = {Bp)\ 
1, 



-)-W. 



(IV.52) 



To find the parameters i?, p,, and W we relate the 
maximums and minimums of the general gray soliton in 
(|IV.52|1 to the maximums and minimums of the trailing 
dip in the DSW of Fig. El Equating maxes and mins 
for the densities, we find 



B 



Po, P 



(IV.53) 



Setting the maximum velocity equal, —W ~ l^fp^ ~ 2, 
we find 



W = 2- 2^. 



(IV.54) 



So the phase speed (|IV.51|) of the gray soliton with the 
parameters (|IV.53|) and (|IV.54p is 



Po, 



equivalent to ljIV.49p . the DSW trailing edge speed. The 
trailing edge of the DSW moves with the phase speed 
of a gray soliton which, in this case, is the sound speed. 
Therefore, similar to the KdV case and, as argued in 
the trailing edge of the DSW can be thought of as a 
modulated train of gray solitons. 

In Fig. El the asymptotic result depicted in Fig. [T7I 
is compared with direct numerical simulation of the NLS 
equation with step initial data showing excellent agree- 
ment. 

The shock profile in Fig. El is valid when 1 < po < 4. 
When Po < 1, the initial data for the Riemann invariants 
r_|. and r_ ljIV.44p is non-decreasing therefore no initial 
data regularization is required and a rarefaction wave so- 
lution exists, see eq. (IIV.29|) and Fig. [201 The rarefaction 
solution for the dispersive and dissipative regularizations 
is the same. 



-W 



D) 

-2- 



B(n-1/n)-W 



FIG. 18: The density and velocity of the gray soliton solu- 
tion 11 V. soil with labeled maximum and minimum values for 
comparison with the trailing dip in the DSW Fig. 1171 



Numerical Solution 



Asymptotic Solution 



5.79 




2.47 



5.79 



2.37 



5.79 



FIG. 19: Numerical solution of the dispersive Riemann prob- 
lem (left) and its asymptotic solution (right) a.t t = 1.5 for 
Po = 2.5 and e — 0.03. The numerically determined trail- 
ing edge speed, calculated from t = 1.2 to t — 1.5, is 1.590, 
approximately equal to y/]}Q = 1.581, the theoretical result 
1IV.49I I. The trailing edge density is 0.175, approximately 
the theoretical value (2 — y/pof = 0.184. 



For the case po > 4, a point of the solution with zero 
density, a vacuum point, is generated [s^l- We label the 
location of the vacuum point with the similarity vari- 
able ^i,. The minimum density of p{x,t) in l)IV.37() is 
Pmin = Ai (since min(— dn(x)2) = — dn(0)2 = —1). Solv- 
ing pmin = Ai = at a vacuum point with the con- 
stants in (|IV.46p (?'!=_ 2, rs = 2, = 4^p^ - 2) 
and using equation (jIV.38p (ri — r2 — ra + r4 = 0) gives 
^3(^1') = — 6. Using the relation (|IV.48|I . the loca- 
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FIG. 20: Numerical solution of the dispersive Riemann prob- 
lem for rarefaction initial data. As e — > 0, the dispersive 
regularization converges strongly to the dissipative regular- 
ization, the rarefaction wave eq. 1IV.29II . The parameters are 
po = 0.2, ito = 2^/po — 2, and t = 0.5. 



tion is 

= 2^-2-2 
my = (Vpo - 1)"^. 



2)A'(m„) 



(IV.55) 



When Po > 4, the characterization of the initial data 
in (|IV.45|I requires that the sign <t in ljIV.37p change at 
the origin (|IV.47p . Since the modulated periodic solu- 
tion l|IV.a7|l and its average PV.40|I depend on a, we 
must determine how a depends on x and t. It is nat- 
ural to assume that the globally conserved momentum 
pu = i|(\E'4'J - **^'^) is smooth. The modulated peri- 
odic wave for the momentum is (recall (|1V.37() ) 

pu = Vp- G\/ A1A2A3. 

Since Y and p are smooth, we consider 

o^/XO^z =^o\^yf^ - 6 - r3(^)| X 
64 

(4Va^ - 2 - r3(e))(4VA^- 6 + r^ii,)). 

The above expression is derived from the rarefaction 
solution H1V.48|I and the relations H1V.38|I . Then, for 
(T\/\\\2\z to be smooth, we require 



Since 



a(0 -sgn(4Vp^-6-r3(C)). 



T-i^iv) = 4Vp^^- 6, 



the sign change occurs exactly at the vacuum point ^ 
^1, ljIV.55|) . Thus, a has the self-similar dependence 



a{x, t) = a{x/t) ^ sgn{x/t - ^i,) 



-1 x/t < 

1 x/t > 

(IV.56) 






FIG. 21: DSW shock profile (solid) and its average (dashed) 
plotted at t = 1 when po > 4 (here po = 10) giving rise to 
a vacuum point at x/t = « 7.4 (see eq. HlV.SSt 'l of zero 
density and infinite velocity. The average velocity V exhibits 
a jump at the vacuum point but the average momentum 4"^ 
does not. This indicates that the variables p and pu are the 
most natural variables for the problem. 



An example DSW profile and its average with a vac- 
uum point at x/t = are shown in Fig. 1211 At the 
vacuum point there is a jump in the average velocity V. 
However, there is no jump in the average momentum tjji/. 

For comparison, Fig. [521 is a plot of the average mo- 
mentum ipi/ (eq. (|IV.40p ) and the modulated DSW 
momentum ^1/ (eq. ljIV.37p ) with the incorrect choice 
a{x,t) = +1. The trailing edge condition 



lim tpi/ = 2po{y/pQ - I) 



(IV.57) 



is not satisfied. Also, the average momentum is not 
smooth at the vacuum point x/t ~ ^y. 

From the plot of the velocity u in Fig. |^ we will 
argue in the next section that vacuum points appear in 
the blast wave experiments considered. 

In summary, the asymptotic solution for a DSW in 
a one-dimensional BEC is a slowly modulated periodic 
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FIG. 22: Average momentum tpv 11 V. 4011 (dashed) and the 
modulated DSW momentum lllV.37l l (sohd) with the rar- 
efaction solution HI V.48II and the incorTect choice of constant 
sign a = +1 and the correct choice a — sgn(^ — ^i,) (dash 
dotted) for po = 10. The boundary condition ljIV.57|l is not 
satisfied and there is a kink at the vacuum point for the in- 
correct a. 



wave expressed as 
p{x,t;e) « Xsix/t) - [Xsix/t) - Xi{x/t)]x 



X dn2[^A3(x/i) - Xi{x/t)e;m{x/t)], 



u{x, t; e) 

m{x/t) 

V{x/t) 
Xiix/t) 

X2{x/t) 

Xsix/t) 



V{x/t) - a{xlt) 



y/Xi{x/t)X2{x/t)X3{x/t) 



p{x,t;e) 



X2{x/t)~Xi{x/t) 
X3{x/t)-Xi{x/t)' 

\[4^-2 + r3{x/t)], 

±[A^-6-r,{x/t)]' 
^[4^+2-r3ix/t)f 
±[4^-2 + r3{x/t)Y 



X - V{x/t)t 



(IV.58) 



where satisfies the implicit equation (|IV.48|I . For 1 < 
Po < 4, a{x/t) = 1. When po > 4, there is a vacuum 
point and a is given in equation (|IV.56|I . 



Here we find that a key difference between dispersive 
and dissipative shock waves is the method of regulariza- 
tion. In the dispersive case, we use initial data regulariza- 
tion. In particular we argue that, with the correct choice 
of initial data, the hyperbolic Whitham equations have a 
smooth solution for all time. Whereas in the dissipative 
case, jump/entropy conditions are employed. Using a 
dispersive regularization, we have determined the behav- 
ior of a fundamental DSW in Bose-Einstein condensates. 
The averaged behavior of the DSW is similar to the clas- 
sical shock case but the speeds and oscillatory behavior 
are different. 



E. Theoretical Explanation of Experiments 

In order to investigate the development of dispersive 
shock waves in 2D and ID, we assume that the conden- 
sate is "prepared" by the 3D evolution, using the re- 
sults of the two experiments with the potentials Vu (eq. 

l|III.2|l ). The state of the condensate 



(jIILT|l ') and Vot (eq. 

at a specific time, t = i (described later), is used as an 
initial condition for a new set of equations in one and two 
dimensions 

2D: «'2D(p,t = 0) = ^{p,z = Q,t = i) 
ID: «'iD(a;,f = 0) = ^(j;, z = 0, t = t), a; >((DV.59) 
*iD(-x,t = 0) = *iD(a;,t = 0). 

Because the anti-trap term in the potentials Vu and Vot 
mainly serves to speed up condensate expansion, we ne- 
glect this term (a^ = 0) when comparing the differing 
dimensional problems. In Fig. 1231 the difference between 
a,- 7^ and = can be seen. Specifically, we solve the 
following three equations numerically in three, two, and 
one dimensions respectively 



3D 



dt 



l£- 



2D 



dt 



2 

1 

2 

f! 

2 



3D 



1 a* 



3D 



3D 



dr 



|*3DP*3D 



2D 



1 (9* 



2D 



dr 



|^2d|'«'20V.6O) 



le- 



ID 



£2 



ID 



dt 



2 Sa 



ID, 



with the initial conditions given in (|1V.59() and 
^3Dir,z,t = 0) = ^{r,z,t = i). 

Comparison of the numerical simulations in different 
dimensions is made by considering the density and 
the appropriate velocity of 5*. In 3D the radial veloc- 
ity is defined to be [arg^'3D(r, 0, t)]r. In 2D, the ve- 
locity is [arg ^'2d(?', i)]r whereas in ID, the velocity is 
[a.rg'i/iY){x,t)]x- The 3D and 2D results are found to be 
barely distinguishable with less than 1% difference be- 
tween them in density, the difference between them can- 
not be seen in Fig. [221 Also, in Figs. 1241 and 1771 it is 
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shown that the three dimensional and one dimensional 
results are qualitatively the same hence the analytical 
studies in ID in section HV Dl arc reasonable approxima- 
tions of the 3D case. 

First we consider the in trap simulation with the po- 
tential Vit (eq. (|III.ip ^. We take the state of the conden- 
sate at time i = 6t {5 ms), just after the laser has been 
applied. Now, the evolution is governed solely by the 
GP equation with an expansion potential. Figure is a 
comparison between the 3D evolution with (a^ = 0.71) 
and without (a^ = 0) the expansion potential. The trail- 
ing edge of the DSW propagates towards the center when 
there is no expansion potential whereas in the full simula- 
tions with the expansion potential, both the leading and 
trailing edges of the DSW propagate radially outward. 
Otherwise, the two simulations are very similar. In order 
to investigate dispersive shock behavior, wc neglect the 
expansion potential, = from now on. We refer to 
this case as the modified in trap experiment. 

The evolution in Fig. [23 shows the development of a 
DSW on the inner ring of high density and another DSW 
on the outer edge of the ring. The density is given on 
the left while the radial velocity is plotted on the right. 
The outer DSW (see t = 2.8 ms) loses its strength and 
vanishes. The DSW shock structure investigated in the 
previous section can clearly be seen in the radial velocity 
plot (compare with the asymptotic solution in Fig. I21|l . 
The appearance of a vacuum point is clear in the velocity 
plot as well. 

In Fig. 1^ the 3D simulation of the modified in trap 
experiment is compared with the corresponding ID sim- 
ulation. Both depict the generation of DSWs and we see 
that they are in good qualitative agreement. The ID den- 
sity grows to about twice the magnitude of the 3D den- 
sity. In the velocity plot, the 3D DSW speeds (trailing 
and leading) are faster than the ID speeds. Therefore, 
the ID analysis performed in section lTVDI is applicable to 
this 3D experiment only in a qualitative sense-it explains 
the basic structure of a DSW in 3D. 

As outlined in section IIVI dissipative and dispersive 
shock waves have quite different properties. In Fig. 1251 
2D simulations of the Euler equations of gas dynamics 
and a BEG are compared. The dissipative regularization 
of the Euler equations was calculated using the finite vol- 
ume package Glawpack ^3 for the conservation laws 

pu 

pt + [pu)r + — = 

{pu)t + {pu"^ + \p^)r + — = 0. 

The same initial conditions for both the dissipative and 
dispersive cases were used. Both simulations depict the 
generation of two shock waves, one on the inner edge 
of the high density ring and another on the outer edge. 
The outer shock vanishes in both cases as it propagates 
into the region of negligible density. Recall that a shock 
wave can only exist when there is a non-zero density on 
both of its sides ^s we have shown, the dissipative 



1*1' (arg*). 




20 30 40 50 60 ( (ms) ^0 30 40 50 60 



FIG. 23: Development and propagation of DSW in the 3D 
in trap experiment with (dashed) and without (solid) the ex- 
pansion potential. The expansion potential does not alter the 
structure of the DSW, just the speeds. The 3D and 2D densi- 
ties are in excellent agreement with a maximum relative error 
of less than 1%. The difference cannot be seen in this graphic. 



1*1' (arg*)r or (arg*)^ 




urn iim 

FIG. 24: Comparison of 3D *3D(r,z = 0, i) (solid) and ID 
'^iD{x,t) (dashed) simulations of DSW in the modified in 
trap experiment showing good qualitative agreement. 
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FIG. 25; Comparison of 2D simulations for BEC (solid) and 
gas dynamics (dashed) . The dissipative regularization used in 
the gas dynamics simulation was calculated using Clawpack 



and dispersive shocks propagate with different speeds. 
However, the densities in the two simulations do agree 
in regions not affected by breaking, the smooth region 
to the right of the DSW. As discussed in the previous 
section, the dissipative and dispersive regularizations are 
the same when there is no breaking. 

For the out of trap experiment, we consider the state 
of the condensate at time t corresponding to about 10 
ms during the experiment or 1 ms after the laser has 
been turned on. This particular i was chosen because 
it represents the time just before breaking occurs in the 
experimental simulations (see Fig. [S]) . As in the previous 
simulations, we neglect the expansion potential (a,, = 0) 
in Vot llTTm . 

In Fig. the development and interaction of 

two DSWs is shown from the evolution of the density 
|*3D(r, 0,t)p and the radial velocity [arg \E'3d(?', 0, i)]r. 
The signatures of two DSWs (Fig. I12|l are easier to see 
in the velocity plots on the right of Fig. I^Hl as the den- 
sity modulations are large. At < = 0.13 ms, two DSWs 
begin to develop on the inner and outer sides of the high 
density ring. At t — 0.33 ms, these two DSWs begin to 
interact giving rise to doubly periodic or multi-phase bc- 



FIG. 26: Three dimensional shock development and interac- 
tion for the modified out of trap experiment. The left plot 
shows the evolution of the density |^'(r, 0,t)|^ and the right 
plot shows the evolution of the radial velocity [arg 'I'(r, 0, t)]r, 
both in the 2 = plane. As with the in trap experiment, the 
maximum absolute difference in density from the 2D evolution 
has a relative error of less than 1%. 



havior. Following this, a DSW front propagates ahead of 
the interaction region. We will examine this interaction 
behavior in more detail in the future. 

In Fig. the 3D and ID simulations are compared 
for the modified out of trap experiment. The agreement 
is quite good over the short time scale considered (1 ms). 
This shows that the initial generation and interaction of 
DSWs in 3D BECs is well explained by the ID approx- 
imation. At later times it is found that the amplitudes 
and speeds diverge roughly similar to what we saw in Fig. 

Figure 1^ shows a comparison of the dissipative regu- 
larization of the 2D Euler equations ljIV.61|) and the dis- 
persive regularization given by the small dispersion limit 
of the Gross-Pitaevskii equation (|IV.60p . The dissipa- 
tive regularization was calculated using Clawpack |l7l |. 
As in the dispersive case, the gas dynamics equations 
develop two shock waves on the inner and outer sides 
of the ring. The inner shock wave propagates as long 
as there is non-zero density on its inner edge. This be- 
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FIG. 27: Comparison of 3D (solid) and ID (dashed) shock 
development. The left plot shows evolution of density; the 
right plot shows the evolution of the velocity. There is good 
qualitative agreement over this short time scale. 
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FIG. 28: Shock evolution in 2D BEG (solid) and the dissipa- 
tive regularization of the 2D Euler equations (dashed). The 
left plot is density and the right plot is radial velocity. The 
dissipative regularization was calculated using Glawpack fl^ . 



havior is predicted by the ID analysis where, with zero 
background density, the solution is always a rarefaction 
(expansion) wave. The structure and speeds of the two 
types of shocks are quite different as expected from the 
analysis given in this work. 



V. CONCLUSION 

We have presented new experimental evidence for dis- 
persive shock waves (DSWs) in a Bose-Einstein conden- 
sate. Numerical simulations of the experiments and com- 
parisons with lower dimensional approximations show 
that the experimentally observed ripples correspond to 
DSWs. A DSW has two associated speeds (front and 
rear of the oscillatory region) and large amplitude oscil- 
latory structure. 

A detailed comparison between classical, dissipative 
shock waves and dispersive shock waves has been un- 
dertaken in this work. On a large scale, a DSW has some 
similarities to a dissipative shock, but there are funda- 
mental and critical differences. Using the notions of av- 
eraging and weak limits, we have shown that a DSW has 
one behavior similar to that of a classical shock-i.e. it 
has a steep front. However, the key difference between 
a DSW and a dissipative shock is in the shock speed 
and its structure. These properties are compared using 
dissipative and dispersive regularizations of conservation 
laws. A dispersive regularization for conservation laws 
gives rise to a weak limit in the sense that one must ap- 
propriately average over the high frequency oscillations 
across a DSW. On the other hand, a dissipative regular- 
ization for conservation laws is a strong limit which, in 
the case of a classical shock wave, converges to a discon- 
tinuity propagating with a speed that satisfies the well 
known jump conditions. 
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